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Mixed-Order  Finite  Elements  for  the  Solution 
of  Three-Dimensional  Electromagnetic  Fields 

Z.J.  Cendes,  D.  Hudak,  D.  Sun 

Electrical  &  Computer  Engineering 
Carnegie-MelJon  University 
Pittsburgh,  PA  15213 

Abstract 

A  new  method  for  modeling  electromagnetic  waves  by  the  finite  element  method  is  presented.  The 
method  is  based  on  a  numerical  formulation  in  which  different  orders  of  polynomials  are  used  to 
approximate  the  three  different  components  of  either  the  electric  or  the  magnetic  field  vectors.  It 
provides  a  reliable  procedure  for  the  finite  element  solution  of  three-dimensional  electromagnetic  field 
problems.  Heretofore,  such  solutions  were  plagued  by  the  presence  of  spurious  modes.  The  new  method 
is  applied  to  the  analysis  of  fields  in  resonant  electromagnetic  cavities. 

Introduction 

The  finite  element  method  is  often  advanced  as  a  useful  numerical  procedure  for  modeling  high- 
frequency  electromagnetic  wave  phenomena.  Applied  at  an  early  date  to  solve  homogeneous  waveguide 
problems,*  the  method  proved  to  be  extremely  accurate  and  reliable  for  these  problems.  However,  when 
the  method  was  applied  to  the  study  of  inhomogeneous  waveguides  and  to  the  solution  of  three- 
dimensional  resonant  cavity  and  scattering  problems,  difficulties  in  the  form  of  ‘spurious*  modes  were 
encountered.  Spurious  modes  are  non-physical  solutions  of  the  electromagnetic  field  equations  that  are 
computed  simultaneously  with  the  correct  physical  solutions. 

Since  the  presence  of  spurious  modes  in  a  numerical  solution  can  destroy  the  validity  of  the  solution, 
much  effort  has  been  directed  at  reducing  or  eliminating  the  unwanted  modes.  The  first  approach, 
originally  suggested  by  Konrad^,  is  to  enforce  the  electromagnetic  field  boundary  conditions  exactly  on 
the  finite  element  approximation  space.  This  procedure  has  been  used  by  Mabaya,  Lagasse  and 
Var.denbulke^  in  the  -  H^  formulation  and  by  Davies,  Fernandes  and  Philippon*  and  by  Rahman  and 
Davies^  in  the  3-component  H  formulation;  all  three  papers  report  only  limited  success  in  eliminating 
spurious  modes  by  this  technique.  Recently,  Koehiba,  Hayata,  and  Suzuki**  have  shown  that  rigorous 
enforcement  of  boundary  conditions  does  indeed  eliminate  spurious  modes  above  the  *air  line*  (i.e.  the 
line  =  1  in  a  ^/k^  versus  k^  plot)  but  does  not  work  in  genera!  below  the  air  line. 

The  second  approach  to  eliminating  spurious  modes  is  to  modify  the  variational  principle  used  to 
approximate  the  fields.  Working  independently,  Winkler  and  Davies^  and  Hara,  Wada,  Fukasawa  and 
Kikuchi*  have  recognized  that  the  spurious  modes  do  not  satisfy  the  zero  divergence  condition  on  the 


electric  or  magnetic  field.  Both  references  suggest  adding  a  penalty  term  proportional  to  the  norm  of  the 
divergence  of  the  field  to  the  governing  variation^  principle.  Unfortunately,  this  procedure  does  not 
eliminate  the  spurious  modes  completely.  However,  as  demonstrated  in  reference  8,  the  spurious  modes 
are  not  stable  with  respect  to  the  amount  of  penalty,  and  can  be  distinguished  from  correct  solutions  by 
plotting  the  finite  element  eigenvalue  spectra  with  respect  to  the  penalty  parameter.  Of  course,  this 
procedure  is  highly  inefficient  and  cumbersome;  each  new  field  problem  must  be  solved  repeatedly  and  all 
of  the  eigenvalues  plotted  in  order  to  identify  the  correct  ones. 

A  different  but  related  procedure  to  reduce  the  number  of  spurious  modes  was  proposed  recently  by 
Konrad.^  References  7  and  8  are  based  on  minimizing  a  functional  derived  from  the  vector  wave  equation 
with  the  addition  of  a  zero-divergence  penalty  term.  Konrad  suggests  using  the  vector  Helmholtz 
equation  instead  and  finds  that  ”...  a  great  number,  though  not  all  of  the  spurious  solutions  are  indeed 
eliminated."  This  result  is  not  surprising  since  the  variational  expressions  derived  for  the  vector 
Helmholtz  equation  and  for  the  vector  wave  equation  with  the  addition  of  a  unit  penalty  term  are 
identical  in  the  case  of  homogeneous  media. 

The  third  approach  to  eliminating  spurious  modes  in  finite  element  solutions  is  to  restrict  the  finite 
element  approximation  functions  to  lie  in  a  reduced  function  space.  In  this  view,  spurious  modes  are  the 
result  of  using  improper  functions  in  the  variational  procedure.  To  ensure  that  only  correct  solutions  are 
generated,  one  must  employ  only  admissible  functions  in  the  finite  element  approximation.  This  is  the 
approach  taken  in  this  report. 

The  use  of  a  restricted  function  space  to  eliminate  spurious  modes  in  finite  element  analysis  was  first 
suggested  by  Hano.^'’  Hano  showed  that,  in  two  dimensional  problems,  spurious  modes  are  completely 
eliminated  by  using  combination  constant-linear  finite  element  ^iproximation  functions.  These  functions 
have  two  interesting  properties:  (1)  their  divergence  is  identically  zero,  and  (2)  they  are  discontinuous  at 
element  boundaries. 

In  this  report,  we  derive  a  set  of  restricted  finite  element  basis  functions  for  the  solution  of  three 
dimensional  electromagnetic  field  problems  and  show  that  only  correct,  physical  solutions  are  obtained 
with  the  new  functions.  As  in  the  case  in  reference  10,  the  functions  reported  here  employ  different  orders 
of  polynomial  approximation  in  different  directions  in  each  element.  However,  in  our  case,  the 
approximation  functions  are  continuous  across  element  boundaries  and  are  not  restricted  to  be  non- 
divergent. 
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Formulation  of  the  EM  Field  EQuationa 
Electromagnetic  wave  propagation  is  governed  by  the  vector  wave  equations 

V  X  —  V  X  E  =  </t*E 

II  **  ® 


(1) 


V  X  ~  V  X  H  =  (2) 

€  ® 
r 

where  and  are  the  relative  permeability  and  relative  permittivity  of  the  material,  respectively,  and 

A  A 

=  u  ft At  the  interface  between  two  dielectrics,  the  tangential  components  of  the  electric  and 
magnetic  fields  must  be  continuous 

X  {E<‘>  -  E^*))  =  0  (3) 


1,  X  =  0 

(4) 

while  the  normal  components  are  discontinuous  as  follows 

o 

It 

V 

1 

• 

(5) 

=  0 

(6) 

In  these  equations,  superscripts  (1)  and  (2)  refer  to  media  1  and  media  2,  respectively. 


In  the  remainder  of  this  report,  we  shall  use  the  electric  field  E  as  the  unknown.  Obviously,  a  similar 
treatment  holds  for  H.  In  terms  of  E,  equation  (4)  becomes 

1  X  (—  V  X  E<‘)  -  —  V  X  E^^h  =  0  (7) 

We  must  therefore  solve  equation  (1)  subject  to  the  interface  conditions  (3),  (5)  and  (7). 

Konrad*  has  shown  that  the  Euler  equation  of  the  functional 

^  ~  \  !  {-  |V  X  E|*  -  (8) 

is  equation  (1)  and  that  the  corresponding  natural  boundary  equation  (7)  is  the  natural  boundary 
condition  for  this  functional.  One  may  therefore  approximate  the  electric  field  E  by  expanding  E  in 
terms  of  basis  functions,  apply  the  essential  boundary  conditions  (3)  and  (5)  explicitly,  and  minimize  F  to 
obtain  E. 

While  the  solution  process  described  above  has  been  widely  reported  in  the  literature^'^^’***'*^,  as  noted 
in  the  introduction,  there  are  serious  problems.  To  eliminate  these  problems,  we  need  to  define  the 
approximation  functions  for  E  more  carefully. 


V  ' 
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A  Basis  for  Curl 

To  solve  for  E  via  equation  (8),  we  must  find  a  legitimate  approximation  for  the  curl  operator.  The 
operator  curl  has  a  domain,  a  nullspace  and  a  range;  it  is  not  enough  —  as  had  been  done  in  the  past  —  to 
approximate  only  the  domain. 

The  nullspace  of  an  operator  is  defined  to  be  the  set  of  functions  that  produce  zero  when  the  operator 
acts  on  it 

N{A)  =  {x:  Ar  =  0  for  all  i}  (9) 

It  is  well  known  that  the  nullspace  of  the  curl  operator  is  provided  by  the  gradient  operator 

N[e\irl)  —  V  4^  (10) 

where  ^  is  an  arbitrary  scalar. 

Let  us  approximate  4  by  finite  basis  functions  a(x,v,z)  over  a  rectangular  parallelepiped.  These  basis 

functions  are  defined  in  Appendix  A  and  result  in  the  expression 

(m,n,p) 

4>(x,y,z)  =  S'  (11) 

^(m,n,p) 

The  polynomial  a  (z,y,z)  is  m’th  order  in  the  x^direction,  n’th  order  in  y,  and  p’th  order  in  z.  Since 
the  derivative  of  an  n’th  order  polynomial  is  (n-l)’st  order,  it  follows  that  the  nullvectors  of  the  curl 
operator  must  have  the  form 

EJ,x,y,x)  =  S' 

E^[x,y,z)  =  S'  (12) 

£;x,v,z)  = 

We  may  write  this  in  the  compact  form 

E  =  tE  (13) 


where 


(m-l,n,p) 


7  =  I  0 


o  0 


(m,n,p-l) 


SJi] 

.  A 


!ti 


The  curl  of  E  is  evaluated  aa 


0 

-d/d  z 

d/dy  E,  ' 

II 

u 

X 

> 

d/d  z 

0 

-d/d  X  E^  ‘ 

1  -d/dy 

Substituting  (13)  into  (16)  gives 

VxE  =  aCE 

where 

1 

d/d  X 

0  E 

!  y(m.n-l,p-l) 

0 

f 

0  ! 

0  = 

0 

0 

0 

_l 

where  the  matrices  D.  are  called  differentiatioD  matrices  and  are  defined  in  Appendix  2. 


We  note  that  the  factorisation  in  equation  (17)  is  not  possible  if  the  same  order  of  polynomial  is  used  to 
approximate  E.  in  all  three  directions. 

The  basis  functions  in  equations  (11),  (13),  and  (17)  provide  a  consistent  representation  of  the  nullspace, 
domain,  and  range  of  the  curl  operator,  respectively.  The  dimension  of  the  nullspace  is  equal  to  the 
number  of  parameters  in  this  is 

Dim{N{eurl))  =  (m+lXn+l)(p+l)  (19) 

The  dimensions  of  the  domain  and  range  spaces  are  equal  to  the  number  of  independent  basis  functions  in 
7  and  0\  these  are 

(20) 

Dim  (Z)(curl))  =  Smnp  +  2mn  +  2fnp  +  2np  +  m  +  n  +  p 

Dim  (R{eurl))  =  3mnp  +  np  +  mp  +mn  (21) 

The  rank  of  the  matrix  C  is  equal  to  the  dimension  of  the  domain  space  minus  the  dimension  of  the 
nullspace 


Rank{C)  =  2mnp  +  mn  +  mp  +  np  —  1 


(22) 
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ComputinK  the  Matrix  Elements 


Substituting  equations  (13)  and  (17)  into  equation  (8)  and  minimizing  with  respect  to  the  coefficients  E 
results  in  the  matrix  equation 

C^KC  E  =  Jfc^M  E  (23) 

where  K  and  M  are  the  matrices 

K  =  /  din  (24) 

M  -  /  dn  (25) 

To  evaluate  the  matrix  elements  in  K  and  M,  it  is  sufficient  to  evaluate  the  integral 

(26) 

Substituting  equation  (A13)  into  (26)  gives 

G  =  (l/I//;  j (3-^(f)  ®  3"(0  (8)  o'iX))^  (^^(f)  (8)  a\()  (8)  o‘(X))djrffdX)  (27) 

Introducing  Kroneckers  identities  (iii)  and  (v)  converts  this  into 

G  =  (1/L/^L,)7<'’)(87<*)(87<‘>  (28) 


7<')  =  y’a(')^(f)  a%)di  (29) 

Thus,  only  the  one  dimensional  T  matrix  is  required  to  evaluate  G.  Numerical  values  of  the  first  two  T 
matrics  are 


T^*)  = 


4  2-1 


—  2  16  2 
30 


-12  4 


To  evaluate  the  finite  element  coefficient  matrix  for  one  element,  one  therefore  needs  to  form  the 
matrices  K  and  M  using  equations  (28)  and  (29),  evaluate  C  by  using  the  values  in  Appendix  B,  and  pre- 
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T 

and  po6t-multiply  K  by  C  and  C,  respectively.  The  contribution  from  each  element  is  combined  with 
the  other  elements  in  the  grid  to  form  a  large,  sparse  matrix  eigenvalu'  problem.  This  eigenvalue  problem 
is  solved  by  using  established  techniques  for  the  eigenvalues  and  eigenvectors  E. 

Computational  Results 

A  computer  program  has  been  developed  based  on  the  above  formulation  to  solve  for  the 
electromagnetic  fields  in  resonant  cavities.  The  program  allows  Tinite  elements  of  mixed  orders  to  be 
assembled  and  solved  for  complex  three-dimensional  geometries.  A  significant  component  of  this  work 
was  the  development  of  a  sparse  matrix  eigenvalue  solution  package  based  on  the  Lanczos  algorithm. 

Some  of  the  computational  results  obtained  from  the  program  are  presented  in  Tables  1-12.  In  these 
tables,  the  approximate  eigenvalues  obtained  by  the  finite  element  method  for  rectangular  parallelpiped 
cavities  are  presented  along  with  the  exact  analytical  values. 

To  begin.  Tables  1  and  2  present  the  eigenvalue  spectrum  for  the  "traditional*  finite  element  solution 
of  the  vector  wave  equation  in  which  each  component  of  the  electric  Held  is  approximated  in  all 
direetiona  by  a  linear  polynomial  (S-node  elements).  Spurious  nodes  are  evident  in  the  spectrum,  both 
below  the  dominant  physical  node  and  above  it.  It  is  obviously  very  difficult  to  distinguish  good  solutions 
from  bad  solutions  with  this  approach. 

Now  consider  the  case  of  constant-bilinear  elements  (m  =  n  —  p  =  1;  4  node  elements).  The 
eigenvalue  spectrums  for  two  different  cavities  obtained  by  using  these  elements  are  presented  in  Tables  3 
and  4.  Notice  that  in  this  case  a  one-to-one  correspondence  exists  between  the  approximate  eigenvalues 
and  the  exact  ones,  with  the  dominant  eigenvalue  being  approximated  reasonably  well.  Although  some 
fairly  large  errors  are  obtained  with  the  high-order  nodes  in  Tables  3  and  4,  these  errors  are  a  result  of 
the  small  number  of  elements  used  in  this  problem.  It  is  significant  that  each  mode  in  these  solutions  can 
be  identified;  greater  accuracy  in  the  eigenvalues  can  be  obtained  when  more  elements  are  employed. 

For  constant-bilinear  elements  the  dimension  of  the  nullspace  isN=2x2x2  =  8.  Corresponding  to 
this  nullity,  eight  lero  eigenvalues  are  found  in  Tables  3  and  4.  Zero  eigenvalues  have  been  computed  in 
this  work  to  confirm  theoretical  predictions;  in  practice,  one  may  save  computer  time  by  computing  the 
positive  eigenvalues  only. 

Solutions  obtained  by  using  linear-biquadratic  elements  (m  =  n  =  p  =  2;  18  node  elements)  are 
presented  in  Tables  5  and  6.  These  solutions  are  seen  to  be  much  mor>>  precise  than  the  constant  bilinear 
element  solutions.  The  reason  for  the  increased  accuracy  is  twofold;  (1)  Higher-order  polynomiab  are 
more  accurate  than  low-order  polynomials,  smd  (2)  linear-biquadratic  elements  are  continuous  across 
element  boundaries  while  the  constant-bilinev  elements  are  not.  A  vector  plot  of  the  field  at  a  height  i 


—  0.35  fcbove  the  bottom  of  the  cavity  is  shown  in  Figure  1.  In  the  case  of  linear-biquadratic  elements, 
the  nullity  of  the  curl  operator  isN>=3x3x3^27.  This  is  equal  to  the  number  of  sero  eigenvalues 


Fisura  It  Plot  of  the  electric  field  for  the  dominafemode  in  a 
cavity  of  dimenrions  ir/S,ir/S,ir/S  at  a  height 
of  X  ■■  0.35  above  the  base 

Finally,  we  report  some  computations  obtained  by  applying  the  flnite  element  method  to  the  vector 
Helmholtt  equation.  The  vector  Helmbolts  equatitm  is  separable  in  a  rectangular  coordinate  system  into 
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three  uncoupled  one-dimensional  Helmholtz  equations.  Not  surprizingly,  the  eigenvalues  of  this  one¬ 
dimensional  equation  are  correctly  computed  as  shown  in  Tables  7  and  8.  However,  if  one  employs  the 
constant-bilinear  4-node  element  in  this  approximation  then  an  incomplete  solution  is  produced.  As 
shown  in  Tables  9-11,  the  eigenvalues  are  correct,  but  only  the  zero-mode  eigenvalues  appear.  Thus, 
solutions  obtained  by  using  constant-bilinear  elements  in  the  Helmholtz  equation  do  not  allow  the 
discontinuities  at  the  element  edges  observed  in  the  solution  of  the  vector  wave  equations. 

When  lineskT-biquadratic  elements  are  used  in  the  vector  Helmholtz  equation,  solutions  in  the  form 
reported  in  Table  12  are  observed.  In  this  ease,  all  eigenvalues  of  the  cavity  are  approximated,  but  the 
zero-mode  eigenvsJues  where  the  solution  is  a  constant  in  one  direction  sje  computed  far  more  accurately 
than  the  non  zero-mode  eigenvalues. 

As  a  final  confirmation  of  the  correctness  of  the  procedure  reported  here,  observe  from  Tables  4  and  7 
that  the  solution  of  the  vector  wave  equation  with  constant-bilinear  elements  and  of  the  vector  Helmholtz 
equation  with  trilinear  elements  are  identical  in  all  modes  to  six  significant  figures.  It  appears  that  these 
two  ways  of  solving  rectangular  cavity  problems  are  numerically  identical.  As  pointed  out  earlier, 
however,  the  advantage  of  solving  the  vector  wave  equation  is  that  it  generalizes  to  modeling  non- 
rectangular  cavaties  and  to  cavities  involving  dielectric  boundaries. 

Conclusions 

Three-dimensional  electromagnetic  field  problems  may  be  solved  by  using  irAed  order  finite  elements  to 
approximate  the  vector  field.  In  this  approach  a  consistent  numerical  approximation  is  made  for  the 
domain,  range,  and  null  spaces  of  the  curl  operator.  Such  a  consistent  approximation  in  the  solution  of 
the  vector  wave  equation  eliminates  the  problem  of  spurious  nodes  that  had  plagued  previous  solution 
procedures. 
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(t)  A(Si{B+C)=A<S)B-^A^C 
(it)  qA  (g)  fiB=Q0{A  (g)  B) 

(Hi)  AB®CD=(A<^C)(B(^D) 

M  (A®B)-^=A-'-^Er'‘ 

(t;)  (A®B)^=A^(^^ 

In  general,  one-dimensional  n’th  order  interpolation  polynomials  are  defined  as 

“'■’.w  =  n  ^ 

where  the  (.  are  the  interpolation  nodes.  Using  the  Kronecker  matrix  product,  two-dimensional 
interpolation  polynomials  are  given  by 

_  3'(")(f){g.3'(’")(0  (AS) 

For  example,  the  second-order  interpolation  polynomial  is  the  nine  element  row  vector 

«i(f)®3(0 

(A6) 

“s(f)"l(0  “3(f)®3(0] 

The  locations  of  the  interpolation  nodes  of  the  polynomiids  are  indicated  in  Figure  2.  During 

computation,  the  nodes  are  stored  most  conveniently  by  columns;  in  this  case,  the  axes  (  and  (  are 
oriented  as  shown  in  Figure  3. 

Now  consider  approximating  an  arbitrary  function  ^((,()  in  terms  of  the  finite  element  approximation 
functions.  We  may  write  this  as 

(A7) 

where  ^  is  an  m  by  n  matrix  of  values  at  the  interpolation  nodes.  Elquation  (A7)  is  converted  into  the 
standard  matrix  form 

(>18) 

by  defining  a  vector  operation  called  ucc  in  the  following  manner; 


VecA  = 


where  A.  is  the  i’th  column  of  the  matrix  A.  The  operator  vec  and  the  Kronecker  product  are  related  by 
the  identity 

vec  ABC  =  (Cf  ®  A)  vecB  (-^10) 


'■f-y 


The  vec  of  a  scalar  is  simply  itself;  therefore,  taking  vec  of  both  sides  of  equation  (A8)  yields 

=  CS^%)  t 

where 

^  =  vec# 


(All) 

{A12) 


Three-dimensional  finite  elements  are  generated  by  an  analogous  procedure  to  that  used  in  two- 
dimensions.  Approximation  functions  for  brick-shaped  elements  are  given  by  the  equation 
3-(n,m.p)  _  (g) a ® 


(A13) 

where  (c>(A)  ve  homogeneous  coordinates  in  the  brick.  Figure  4  presents  the  node  numbering  scheme  for 
the  second  order  case. 
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APPENDIX  -  B 
Differential  Matrices 

One-dimensional  differentiation  matrices  are  defined  by  the  equation 


(Bl) 


where  is  the  n  by  n+1  differentiation  matrix.  The  polynomials  *\f)  in  equation  (Bl)  are  of 
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CURL-CURL-trilinear  functions (8-node  elements) 

Lx=Ly=Lz=pi/3.  (pl=3. 1415926535) 

number  of  nodes  in  x-direction  =  12 

"  y-direction  =  12 

"  z-direction  =  12 

hx=hy=hz=pi/9. 


approx-eigenvalue  modes  actual-eigenvalue 


1.96968 

9.08752 

9.08752 

9.08752 

19.69684 

19.69684 

19.69684 

19.69684 

19.69684 

19.69684 

28.25468 

28.25468 

28.25468 

28.56042 

28.56042 

49.24210 

59.09051 

59.09051 

59.09051 

59.09051 

59.09051 

59.09051 

59.09051 

59.09051 

59.09051 

59.09051 

59.09051 

59.09051 

60.83625 

60.83625 

60.83625 

67.95409 

67.95409 

67.95409 
83.71156 
83.71156 
83.71156 
98.48419 
98.48419 
98.48419 
98.48419 
98.48419 
98.48419 


(0,1,1)  18.00000 
(1,0,1) 

(1,1,0) 


(0,1,2)  45.00000 

(0.2,1) 

(1.0,2) 

(2,0,1) 

(1.2,0) 

(2,1,0) 


(0,2,2)  72.00000 

(2,0,2) 

(2,2,0) 


Table  1.  Eigenvalues  of  the  vector  wave  equation  (curl-curl  operator)  in  a  rectangular 
cavity  obtained  by  using  trilinear  approximation  functions.  Note  that  many  of  the 
approximate  eigenvalues  computed  do  not  correspond  to  any  actual  physical  mode. 
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CURL-CURL-trilinear  functions (S-ncxSe  elements) 

Lx=Lz=pi/3.  ,  Ly=pi/4.  (pi=3. 1415926535) 

number  of  nodes  in  x-direction  =  12 

"  y-direction  =  12 

"  z-direction  =  12 

hx=hz=pi/9.  ,  hy=pi/12. 


;  approx-eigenvalue 

2.37084 

I  9,83625 

11.72376 
11.72376 
19.69684 
19.69684 
27.35672 

\  27.35672 

27.35672 
27.35672 
30.68483 

34.39606 

34.39606 
35.81914 
36.22030 
58.39643 
59.09051 

^  .  59.09051 

59.09051 

I  59.09051 

66.12458 
66.12458 
66.75040 
66.75040 
66.75040 

I  66.75040 

75.34928 
75.34928 
91.37144 

97.38992 

97.38992 

37 . 38992 

97.38992 
98.38694 
98.48419 
98.48419 


Table  2.  Similar  results  to 


i 


modes  actual-eigenvalue 


(1,0,1) 

18.00000 

(0,1,1) 

25.00000 

(1,1,0) 

(1,0,2)  45.00000 

(2,0,1) 


(0,1,2)  52.00000 

(2,1,0) 


(0,2,1)  73.00000 

(1,2,0) 

(2,0,2)  72.00000 

Table  1  but  with  a  different  cavity  shape. 


CURL-CURij-bilinear  elements  (4-ncxle  elements) 

Lx=Lz=pi/3.  ,  Ly=pi/4.  (pi=3. 1415926535) 
number  of  nodes  in  x-direction  =  11 

"  y-direction  =  11 

"  z-direction  =  ll 

hx=hz=pi/9.  ,  hy=pi/12. 

approx-eigenvalue  modes  actual-eigenvalue 


0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

19.69684 

(1,0,1) 

18.00000 

27.35672 

(0,1,1) 

25.00000 

27.35672 

(1,1,0) 

37.20514 

37.20514 

(1,1,1) 

34.00000 

59.09052 

(1.0,2) 

45.00000 

59.09052 

(2,0,1) 

66.75040 

(0,1,2) 

52.00000 

66.75040 

(0,2,1) 

76.59882 

76.59882 

(1,1,2) 

61.00000 

76.59882 

76.59882 

(2,1,1) 

97.38992 

(0,2,1) 

73,00000 

97.38992 

(1,2,0) 

98.48419 

(2,0,2) 

72.00000 

107.23834 

107.23834 

(1,2,1) 

82.00000 

115.99249 

115.99249 

(2,1,2) 

88.00000 

Table  3.  The  eigenvalues  of  the  vector  wave  equation  obtained  by  using  constant- 
bilinear  approximation  functions.  The  zero  modes  may  be  ignored.  Notice  that  all 
actual  physical,  eigenvalues  are  correctly  approximated  by  this  procedure  although 
errors  in  the  values  are  large  for  the  higher-order  modes. 


CURL-CURL-bilinear  functions {4-ncKSe  elements) 


Lx=Ly=Lz=pi/3.  (pi=3. 1415926535) 
number  of  nodes  in  x-direction  =  11 

"  y-direction  =  11 

"  z-direction  =  11 

hx=hy=hz=pi/9. 


approx-eigenvalue 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

o.oooco 

0.00000 

modes 

actual-eigenvalue 

19.69684 

(0,1,1) 

18.00000 

19.69684 

(1.0,1) 

19.69684 

(1,1,0) 

29.54526 

29.54526 

(1,1,1) 

27.00000 

59.09052 

(0,1,2) 

45.00000 

59.09052 

(0.2,1) 

59.09052 

(1,0,2) 

59.09052 

(2,0,1) 

59.09052 

(1,2,0) 

59.09052 

(2,1,0) 

68.93893 

68.93893 

(1,1,2) 

54.00000 

68.93893 

68.93893 

(1,2,1) 

68.93893 

68.93893 

(2,1,1) 

98.48419 

(0.2,2) 

72.00000 

98.48419 

(2,0,2) 

98.48419 

(2,2,0) 

108.33261 

108.33261 

(1,2,2) 

81.00000 

108.33261 

(2,1,2) 

108.33261 

108.33261 

108.33261 

(2,2,1) 

Table  4.  Similar  results  to  Table  1  but  with  a  different  cavity  shape 


CURL-CURL-linear-biquadratic (18-no(3e  elements) 


Lx=Ly=Lz=pi/3.  (pi=3. 1415926535) 

number  of  nodes  in  x-direction  =  18 

"  y-direction  =  18 

"  z-direction  =  18 

nx=ny=hz=pi/9 . 


approx-eigenvalue  modes  actual-eigenvalue 

0.00000  -  - 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 

0.00000 


18.02846 

(0,1,1) 

18.00000 

18.02846 

(1,0,1) 

18.02846 

(1,1,0) 

27.02822 

(1,1,1) 

27.00000 

27.02822 

45.80279 

(0,1,2) 

45.00000 

45.80279 

(0,2,1) 

45.80279 

(1,0,2) 

45.80279 

(2,0,1) 

45.80279 

(1,2,0) 

45.80279 

(2,1,0) 

54.06912 

(1.1,2) 

54.00000 

54.06912 

54.06912 

(1,2,1) 

54.78579 

54.78579 

(2,1,1) 

54.78579 

73.57712 

73.57712 

73.57712 

78.67638 

78.67638 


(0,2,2) 

(2,0,2) 

(2,2,0) 


72.00000 


CURL-CURL-linear-biquadratic (18-node  elements) 


Lx=pi/2.  ,  Ly=pi/3.  ,  LZ=pi/4.  (pi=3. 1415926535) 
number  of  nodes  in  x-direction  =  18 

"  y-direction  =  18 

"  z-direction  =  18 

hx=pi/6.  ,  hy=pi/9.  ,  hz=pi/16. 


approx-eigenvalue 

0.00000 
0.00000 
o.ooooc 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
0.00000 
13.02056 
20.03163 
25.03953 
.  25.36470 
29.02859 
29.03043 
32.37577 
40.00222 
40.79489 
41.36763 
45.48986 
52.50093 
52.81386 
53.13903 
55.97965 
56.78009 
61.47941 
61.49460 
62.74550 
68.25079 
69.40821 
73.26419 


modes  actual-eigenvalue 


(0,1,1) 

(1,0,1) 

(0,1,1) 

(2,1,0) 

(1,1.1) 

(2,0,1) 


(1,2,0) 

(2,1,1) 

(3,1,0) 

(0,2,1) 

(2,2,0) 

(3,0,1) 


(1,2,1) 

(3,1,1) 


13.00000 

20.00000 

25.00000 

29.00000 

32.00000 


40.00000 

41.00000 

45.00000 

52.00000 


56.00000 

61.00000 


(2,2,1)  68.00000 
(1,0,2) 

(3,2,0) 


72.00000 


LAPLACIAN-trilinear  functions  (8-no<3e  elements) 


Lx=Ly=Lz=pi/3.  (pi=3. 1415926535) 
number  of  nodes  in  x-direction  =  4 

"  y-direction  =  4 

"  z-direction  =  4 


hx=hy=hz=pi/9. 

approx-eigenvalue 

(19.69684,0.0000000)  2 

(29.54527,0.0000000)  3 

(59.09052,0.0000000)  5 

(59.09051,0.0000000) 
(68.93894,0.0000000)  6 

(68.93893,0.0000000) 
(68.93893,0.0000000) 
(98.48419,0.0000000)  8 

(108.3326,0.0000000)  9 

(108.3326,0.0000000) 
(108.3326,0.0000000) 


modes  actual-eigenvalue 

(0,1,1)  (18.00000,0.0000000) 

(1.1.1)  (27.00000,0.0000000) 

(0,1,2)  (45.00000,0.0000000) 

(0,2,1) 

(1.1.2)  (54.00000,0.0000000) 

(1,2,1) 

(2,1,1) 

(0,2,2)  (72.00000,0.0000000) 

(1.2.2)  (81.00000,0.0000000) 

(2,1,2) 

(2,2,1) 


Table  7.  Eigenvalues  of  the  vector  Helmholtz  equation  (Laplacian  operator) 
obtained  by  using  trilinear  approximation  functions.  Note  that  all  actual 
physical  eigenvalues  are  correctly  approximated  by  this  procedure. 


LAPLACIAN-trilinear  functions (S-node  elements) 


Lx=L2=pi/3.  ,  Ly=pi/4.  (pi=3 .1415926535) 
number  of  nodes  in  x-direction  =  4 
"  y-direction  =  4 

"  z-direction  =  4 

hx=hz=pi/9.  ,  hy=pi/12. 

approx-eigenvalue  modes  actual-eigenvalue 

(27.35672,0.0000000)  (0,1,1)  (25.00000,0.0000000) 

(37.20514,0.0000000)  (1,1,1)  (34.00000,0.0000000) 

(66.75041,0.0000000)  (0,1,2)  (52.00000,0.0000000) 

(76.59882,0.0000000)  (1,1,2)  (61.00000,0.0000000) 

(76.59882,0.0000000)  (2,1,1) 

(97.38994,0.0000000)  (0,2,1)  (73.00000,0.0000000) 

(107.2384,0.0000000)  (1,2,1)  (82.00000,0.0000000) 

(115.9925,0.0000000)  (2,1,2)  (88.00000,0.0000000) 

(125.8409,0.0000000)  (0,1,3)  (97.00000,0.0000000) 

(136.7836,0.0000000)  (0,2,2)  (100.0000,0.0000000) 

(146.6320,0.0000000)  (2,2,1)  (109.0000,0.0000000) 

(146.6320,0.0000000)  (1,2,2) 


Table  8.  Similar  results  to  Table  7  but  with  a  different  cavity  shape 
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LAPLACIAN-bilinear  functions (4-no<ae  elements) 

Lx=Ly=Lz=pi/3.  (pi=3. 1415926535) 

number  of  nodes  in  x-direction  =  4 

"  y-direction  =  4 

"  z-direction  =  3 

hx=hy=hz=pi/9. 


approx-eigenvalue 

modes 

actual-eigenvalue 

(19.69684,0.0000000) 

( 19 . 69684 , 0 . 0000000 ) 

( 19 . 69684 , 0 . 0000000 ) 

(0,1,1) 

(18 . 00000 , 0 . 0000000 ) 

(59.09051,0.0000000) 

(59.09051,0.0000000) 

(59.09051,0.0000000) 

(0,1,2) 

(45.00000,0.0000000) 

(59.09051,0.0000000) 

(59.09051,0.0000000) 

(59.09051,0.0000000) 

(0,2,1) 

(45.00000,0.0000000) 

(98.48419,0.0000000) 

(98 . 48419 , 0 . 0000000 ) 
(98.48419,0.0000000) 

(0,2,2) 

(72.00000,0.0000000) 

Table  9.  Eigenvalues  of  the  vector  Helmholtz  equation  obtained  by  using  constant- 
bilinear  approximation  functions.  Notice  that  only  the  0-mode  eigenvalues  are 
produced. 
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LAPLAC I AN-bi linear  functions  (4-no<le  elements) 

Lx=Lz=pi/3.  ,  Ly=pi/4.  (pi=3. 1415926535) 
number  of  nodes  in  x-direction  =  4 

"  y-direction  =  4 

"  z-direction  =  3 

Iix=hz=pi/9.  ,  hy=pi/12. 


approx-eigenvalue 

modes 

actual-eigenvalue 

(27.35672,0.0000000) 

(27 . 35672 , 0 . 0000000) 
(27.35672,0.0000000) 

(0,1,1) 

(25 . 00000 , 0 . 0000000 ) 

(66.75039,0.0000000) 

(66.75039,0.0000000) 

(66.75039,0.0000000) 

(0,1,2) 

( 52 . 00000 , 0 . 0000000) 

(97.38992,0.0000000) 

(97.38992,0.0000000) 

(97.38992,0.0000000) 

(0,2,1) 

(73.00000,0.0000000) 

(136.7836,0.0000000) 

(136.7836,0.0000000) 

(136.7836,0.0000000) 

(0,2,2) 

( 100 . 0000 ,0.0000000) 

Table  10.  Similar  results  to  Table  9  but  with  a  different  cavity  shape 


LAPLACIAN-bilinear  functions (4-node  elements) 


Lx=Lz=pi/3.  ,  Ly=pi/4.  (pi=3. 1415926535) 
number  of  nodes  in  x-direction  =  5 

"  y-direction  =  7 

"  z-direction  =  6 

hx=hz=pi/15.  ,  hy=pi/24. 

approx-eigenvalue  modes  actual-eigenvalue 

(25.66859,0.0000000)  (0,1,1)  (25.00000,0.0000000) 

(25.66859,0.0000000) 

( 25 . 66859 , 0 . 0000000) 

(25.66859,0.0000000) 

(25.66859,0.0000000) 

(57.30186,0.0000000)  (0,1,2)  (52.00000,0.0000000) 

^7. 30186, 0.0000000) 

(57.30186,0.0000000) 

(57.30186,0.0000000) 

(57.30186,0.0000000) 

(79.33304,0.0000000)  (0,2,1)  (73.00000,0.0000000) 

(79.33304,0.0000000) 

(79.33304,0.0000000) 

(79.33304,0.0000000) 

(79.33304,0.0000000) 

(110.9663,0.0000000)  (0,2,2)  (100.0000,0.0000000) 

(110.9663,0.0000000) 

(110.9663,0.0000000) 

(110.9663,0.0000000) 

(110.9663,0.0000000) 


Table  11.  Similar  results  to  Table  10  but  with  a  finer  discretization. 


LAPLACIAN-linear-biquadratic  functions (IB-node  elements) 

Lx=Ly=L2=pi/3.  (pi=3. 1415926535) 
number  of  nodes  in  x-direction  =  4 

"  y-direction  =  7 

"  2-direction  =  7 

hx=hy=h2=pi/9. 


approx-eigenvalue 

modes 

actual-eigenvalue 

18.02849 

(0,1,1) 

18.00000 

27.87692 

(1,1,1) 

27.00000 

45.80283 

(0,1,2) 

45.00000 

45.80284 

(0,2,1) 

55.65121 

(1,1,2) 

54.00000 

55.65127 

(1,2,1) 

67.27062 

—  ?  — 

73.57719 

(0,2,2) 

72.00000 

83.42559 

(1,2,2) 

81.00000 

91.08446 

(0,3,1) 

90.00000 

91.08447 

(0,1,3) 

95.04494 

—  7  — 

95.04499 

—  7  — 

Table  12.  Eigenvalues  of  the  vector  Helmholtz  equation  obtained  by  using  linear 
biquadratic  approximation  functions.  In  this  case,  all  eigenvalues  are  approxi 
mated,  but  the  0-mode  eigenvalues  are  far  more  accurate  than  the  non  0-mode 
eigenvalues. 


LAPLACIAN-linear-biquadratic (18-node  elements) 


Lx=pi/2.  ,  Ly=pi/3 

.  ,  Lz=pi/4. 

(pi=3. 1415926535) 

number  of  nodes  in 

x-direction  = 

4 

fi 

y-direction  = 

)»7 

n 

z-directlon  = 

hx=pi/6.  ,  hy=pi/9 

.  ,  hz=pi/12. 

approx-eigenvalue 

modes 

actual-eigenval ue 

25,03956 

(0,1,1) 

25.00000 

29.41665 

(1,1,1) 

29.00000 

46.92494 

(2,1,1) 

41.00000 

52.81389 

(0,2,1) 

52.00000 

57.19098 

(1,2,1) 

56.00000 

68.81034 

(0,1,2) 

68.00000 

74.41617 

(1,1,2) 

72.00000 

74.69928 

(3,1,1) 

61.00000 

78.79325 

(2,2,1) 

68.00000 

96.30156 

(2,1,2) 

84.00000 

96.58467 

(3,2,1) 

88.00000 

98.09551 

(0,3,1) 

97.00000 

102.1905 

(0,2,2) 

100.0000 

Table  13.  Similar  results  to  Table  12  but  with  a  different  cavity  shape 
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